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We test the finite density algorithm in the canonical ensemble which combines the HMC update with the 
accept /reject step according to the ratio of the fermion number projected determinant to the unprojected one as 
a way of avoiding the determinant fluctuation problem. We report our preliminary results on the Polyakov loop 
in different baryon number sectors which exhibit deconfinement transitions on small lattices. The largest density 
we obtain around T c is an order of magnitude larger than that of nuclear matter. From the conserved vector 
current, we calculate the quark number and verify that the mixing of different baryon sectors is small. 



1. Motivation 

The phase structure of QCD is relevant in the 
description of various phenomena: from subtle 
modifications of the cross-sections in high energy 
collisions of nuclei to exotic states of nuclear mat- 
ter in neutron stars. Even though we have a good 
understanding of how to perform simulations at 
zero baryon density and finite temperature, the 
simulation with non-zero chemical potential has 
been problematic. 

To deal with the complex phase that the chem- 
ical potential introduces to the fermion determi- 
nant, a number of techniques have been proposed: 
various reweighting methods, imaginary chemi- 
cal potential, etc. They allow us to explore the 
phase structure at small values of the chemical 
potential and temperatures close to the critical 
temperature. All these methods have as a start- 
ing point the grand canonical partition function. 
We will present here an approach based on the 
canonical partition function which alleviates the 
sign problem, the overlap problem and the deter- 
minant fluctuation problems |1I2| . Some prelimi- 
nary results will also be reported. 

2. Canonical partition function 

The QCD canonical partition function with a 
fixed quark number can be constructed from the 
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projection of the quark number from the deter- 
minant |3I2| . The probability of having k quark 
loops more than antiquark loops wrapping around 
the time direction in the loop expansion of a de- 
terminant is represented by the Fourier transform 

1 r 27r 

p k = 7T d0e^ fc detM . (1) 
27T Jo 

where det is the fermion determinant of the 
quark matrix with an additional U(l) phase 
e %4> I e -i4> j n ftie time forward/backward links on 
one time slice. Therefore, the canonical partition 
function can be written as 

Z c (k,T) = J DU e~ SG ^e~^ k detM . (2) 

One can also derive this from the Fourier trans- 
form of the grand canonical partition function 4 

Z c (k,T) = J 2 J ^e-** k Zac(ji 9 = ^T,T) 

To see this we can use the fugacity expansion of 
the grand canonical partition function 

12V 

Z GC ( N ,T)= J2 Z c (k,T)e^ T (3) 

k=-12V 

where k = n q — riq is the quark number and [i q 
is the quark chemical potential. Using the usual 
expression for Lattice QCD grand canonical par- 
tition function, one obtains the same expression 
as in Eq. ©. 
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For Wilson fermions we can write M^U) = 
M(U$) = M(fi = icj)T, U) in terms of the stan- 
dard fermion matrix 

[M(U)] m ,n = 8 m . n - k(1 - 7 M )[/ M (m)(5 m+M ,„ 
- «(l + 7 Al )lJt(n)(5 



and some modified gauge fields 



_ f £^(n)e-<* n 4 = 7V t , M = 4 
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(4) 



This should not be viewed as a change of the field 
variables but rather as a convenient way of writ- 
ing the partition function. We should also note 
that although originally the phase e~^ T appears 
on all timclikc links, we can accumulate it on the 
last time slice through a change of variables. 

In order to evaluate the canonical partition 
function in Eq. @ we need to replace the contin- 
uous Fourier transform with a discrete one. Using 
the discrete Fourier transform of the determinant 

N-l 

det k M{U) = ^Y, det M *n( u ) 



n=0 



with <f) n = ^jp- , we write the partition function as 



Z c (k) = J BUe- SG{u) det k M(U) 



(5) 



This is the partition function that we will try to 
evaluate. The parameter TV defines the Fourier 
transform. In the limit N — > oo we recover the 
original partition function. For finite N the parti- 
tion function will only be an approximation of the 
canonical partition function. Using the fugacity 
expansion we can show that 

oo 

Z c {k)= Z c (k + mN) 

m— — oo 

If Fb denotes the baryon free energy we expect 
that 

Z C (N B ) oc e -l"B|iWr 

in the confined phase. If this holds true then 
Z c (k) » Z c (k) for k < N. However, this is 
an expectation that needs to be checked in our 
simulations. 



3. Algorithm 

To evaluate Zc{k) by Monte-Carlo techniques 
we need the integrand in J^J) to be real and posi- 
tive. Using the 75 hermiticity of the Wilson ma- 
trix it is easy to prove that det Af (J/^,) is real. 
However, the integrand in J^} is the Fourier trans- 
form of the Wilson determinant. For detfcAf 
to be real we need that detM^ = detM_^. 
From charge conjugation symmetry we know that 
(det A/0 ) = (det M_^). Unfortunately this rela- 
tion doesn't hold configuration by configuration. 
We are thus forced to remove the complex phase 
from detfcAf. 

To avoid dealing with the complex phase we 
will generate an ensemble using |Re det^Af | as our 
probability function and we will then reintroduce 
the complex phase in observables. This is the 
standard way to deal with a complex phase in 
the weight function. It is also the source of the 
infamous sign problem. This might limit us in 
fully exploring the parameter space. 

The focus of the algorithm will then be to gen- 
erate an ensemble with the weight 



W(U) 



S G (U) 



N-l 

— Y cos(0„fc)detAf „(Z7) 
n=0 



It was pointed out in [2] that taking the Fourier 
transform after the Monte Carlo simulations with 
different <p n leads to an overlap problem. To avoid 
this problem, it was stressed [21 1| that it is es- 
sential to perform the Fourier transform first to 
lock into one particular quark (or baryon) sector 
before the accept/reject step. The general ap- 
proach is to break the updating process into two 
steps. In this case, one first uses an update for an 
approximate weight W'(U), and then use an ac- 
cept/reject step to correct for the approximation. 
The efficiency of this strategy depends on how 
good the approximation is. For a poor approxi- 
mation the acceptance rate in the second step will 
be very small. 

One possible solution would be to update the 
gauge links using a heatbath for pure gauge ac- 
tion. In that case W'(U) = e~ Sa ( u K However, it 
is known that such an updating strategy is very 
inefficient since the fermionic part is completely 
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disregarded in the proposal step and the deter- 
minant, being an extensive quantity, can fluctu- 
ate wildly from one configuration to the next in 
the pure gauge updating process To avoid 
the overlap problem and the leading fluctuation 
in the determinant, it is proposed pQ to update 
using the HMC algorithm in the first step. For 
simplicity we will be simulating two degenerate 
flavors of quarks. This allows us to use the stan- 
dard HMC update. In this case we have 

W'(U) = e ~ SG(c/) det M(U) 

In the second step we have to accept/reject 
with the probability 

P acc = min{l,u(U')/w(U)} 

where u> is the ratio of the weights 

W(U) |Redet fc M[J7]| 



,{U) = 



W'{U) 



det M[U] 



Since detkM[U] and detM[J7] are calculated on 
the same gauge configuration, we can write the 
determinant ratio as 



u(U) = — 



N-l 



cos(0„/fc)e T ''( logM ^- logM ) 



n=0 



(6) 



The leading fluctuation in the determinant from 
one gauge configuration to the next is removed 
by the Tr log difference of the quark matrices of 
M^, n ({7) and M(U). This should greatly improve 
the acceptance rate compared to the case with 
direct accept/reject based on the |RedetfcM([/)| 
itself. 

4. Triality 

The canonical partition function (J2J) has a Z% 
symmetry that is not present in the grand canon- 
ical partition function 7.. Under a transforma- 
tion U — > Utj, where U$ is defined as in @J and 
4> = ±^ we have 



det k M[U] -> det k M[U ±2 ^] 



det k M[U). 



We see that when A: is a multiple of 3 this transfor- 
mation is a symmetry of the action. Incidentally, 
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Figure 1. Polyakov loop argument as a function 
of simulation time. 



if this symmetry is not spontaneously broken the 
relation above guarantees that the canonical par- 
tition function will vanish for any k that is not a 
multiple of 3. 

For Zc, our approximation of the canonical 
partition function, to exhibit this symmetry we 
need to choose the parameter N that defines the 
discrete Fourier transform to be a multiple of 3. 

If N is a multiple of 3 we have 

W{U)=W{U ±¥ ). 

However, the HMC weight W does not respect 
this symmetry since det M is not symmetric un- 
der this transformation. In this case, our algo- 
rithm can become frozen for long periods of time. 
In Fig. 2]we show the argument of the Polyakov 
loop as a function of the simulation time. Under 
a 2T3 transformation U — > U±m. the argument of 

the Polyakov loop gets shifted by The HMC 

update strongly prefers the sector and when a 
tunneling occurs to one of the other sectors the 
accept/reject step will reject any change for a long 
period of time. Consequently the algorithm will 
be very poor in sampling the other sectors. 

To alleviate this problem we will introduce a 
2T3 hopping Since the weight W is symmet- 
ric under this change we can intermix the regu- 
lar updates with a change in the field variables 
U — > E/±2^, where we will choose the sign ran- 
domly (equal probability for each sign) to satisfy 



4 



Table 1 Table 2 

The simulation parameters for our runs. Average of real phase factor, (Rea)w- 

(3 I q(fm) I m^MeV) I y^fm" 3 ) I r(McV)~| I I k=0 j k=3 | k=6 ~ 

5.1 0.35(3) 870(90) 0.36(12) 140(14) /3=5.1 1.00(0) 0.48(6) 0.17(10) 

5.2 0.31(2) 920(90) 0.52(17) 160(16) /3=5.2 1.00(0) 0.61(8) 0.82(5) 

5.3 0.24(2) 1050(100) 1.1(3) 205(20) g=5j 1.00(0) 0.98(1) 0.96(3) 



the detailed balance. This will decrease our ac- 
ceptance by a factor of roughly 3 but will sample 
all the Z% sectors in the same manner. 

5. Simulation details 

In order to simulate the ensemble given by the 
weight W(U) we need to evaluate the Fourier 
transform of the determinant. Since the calcula- 
tion of the determinant is quite time consuming 
we will need to use N as small as possible. There 
is a proposal that employs an estimator for the 
determinant [§] but in this preliminary study we 
choose to calculate the determinant exactly us- 
ing LU decomposition. We were thus limited to 
running on very small lattices, i.e. 4 4 , and with 
N = 12. For each (3 we run three simulations: 
one for k — 0, one for k = 3 and one for k = 6. 
They correspond to 0, 1 and 2 baryons in the box. 

In order to get a reasonable box size we had 
to simulate at large lattice spacings. We used 
(3 = 5.1, 5.2 and 5.3 with fixed k = 0.158. 
The relevant parameters can be found in Table 
□ The lattice spacing is determined using stan- 
dard dynamical action on a 12 4 lattice for the 
same values of (3 and k. For the HMC part of 
the update we used trajectories of length 0.5 with 
At = 0.01. We adjusted the number of such tra- 
jectories between two consecutive accept/reject 
steps so that the acceptance stays within the 15% 
to 30% range. The configurations were saved af- 
ter 10 accept/reject steps. We collected about 
100 such configurations for each run. 

6. Observables 

As we mentioned earlier we need to reintroduce 
the phase when we compute the observables. For 



an observable O we have the following expression 



where the subscript Zc stands for integration 
over the canonical ensemble given by J5J, W 
stands for the average over the weight W and 

det fe M 
a = |Redet fe M| 

is the reweighting phase. It is easy to see that the 
phase a should average, over an infinite ensemble, 
to a real number. This is direct consequence of 
the charge conjugation symmetry. The real part 
of this phase will always be ±1. This gives us a 
simple criterion to decide whether we have a sign 
problem; all we need to do is to check whether the 
phase has predominantly one sign. If the phase 
comes with plus and minus signs with close to 
equal probability then we have a sign problem. 
For our runs the data is collected in Tabled We 
see that for (3 = 5.1 and k = 6 we are approaching 
the sign problem region. 

An interesting observable to measure is the 
Polyakov loop. Although the average value is 
expected to be zero due to the Z3 symmetry, 
the distribution in the complex plane can tell us 
whether or not the symmetry is spontaneously 
broken. Usually the spontaneous breaking of the 
Z3 symmetry is associated with deconfinement. 
In Fig. [5]we show the histograms of the Polyakov 
loop for our runs. We see that as we go from 
(3 = 5.2 to P = 5.3 the symmetry is broken even 
in the k = sector. This is the usual finite tem- 
perature phase transition at zero baryon density. 
According to TableHthis happens somewhere be- 
tween 160 MeV and 205 MeV. This is somewhat 
higher than expected for dynamical simulations 
but can be explained in terms of our large quark 
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Figure 2. Polyakov loop distribution. 



masses. More interestingly, we see that under the 
critical temperature the transition still occurs as 
we increase the number of baryons. It is also ap- 
parent that as we lower /3 (/? = 5.1), we are ap- 
proaching a region where even the Ng — 1 sector 
is confined. 

Assuming that the result of Nb = l(fc = 3) at 
(3 = 5.1 is close to the phase transition line, the 
phase transition for T = 140 MeV occurs at the 
baryon density of 0.36 fm -3 which is about 2.3 
times the nuclear matter density. We see from 
Table 1 that there does not seem to be a problem 
of reaching a density an order of magnitude higher 
than that of the nuclear matter around T c . For 
example, for Nb — 2(fc = 6) at j3 = 5.3, the 
baryon density is 2.2(6) furi 3 which is 14 times 
the nuclear matter density. However, we should 
not take these numbers seriously in view of the 
small volume and large quark mass in the present 
calculation (since the pion mass is 870 MeV our 
quark mass is close to the strange quark one). 

While it is nice to see that our expectation for 
the phase transition are confirmed in the gluonic 
sector, one would like to explore this transition in 
the fermionic sector too. It is important to note 
that the fermionic observables have to be evalu- 
ated in a non-standard manner. For example, if 
we are to study the fermionic bilinears we get 



1 1 
~Zc~N 



where 



Tr^M" 1 = 1 ^e-^-TiTA^- 1 . 

We see that for fermionic variables, we need to 
do a Fourier transform on the observable too. 
We have measured the chiral condensate in our 
runs in the hope that it will signal chiral sym- 
metry restoration. Unfortunately, it seems that 
the quarks are too heavy for us to get any signal. 
For example at (3 = 5.1 and k — 0.158 we got 
(■ipip) = -3.612(1), -3.598(5) and -3.601(2) for 
Nb = 0,1 and 2 respectively. We suspect that we 
need to lower the quark mass to be able to see any 
signal in the chiral condensate. Also, since we are 
working with Wilson fermions, the signal might 
be obscured by the mixing between tpif) and unity 
and other lattice artifacts especially at these large 
lattice spacings. 

7. Sector mixing 

It is important to note that all the calculations 
carried out above are using an approximation of 
the canonical partition function. As noted before 
it is an assumption on our part that Zc ~ Zc- 
To check this we measured the conserved charge 
for Wilson fermions 

Q(t) = -k^2 $(s)E/k(a;)(l- 74)^(0; + *) 

X 

-^ + £)C/ 4 t (x)(l + 74)V(^)]. 

The idea is that if we were to compute the con- 
served charge using the exact canonical ensem- 
ble Zc we would measure exactly the number of 
quarks we put it (Q(t)) z a (k) — However, in 
the discrete case we have 



(Q(*)> 



Zc(k) 



J2 m {k + mN)Z c (k + mN) 
J2 m Z c (k + mN) 



This is a weighted average over the various sec- 
tors allowed by our discretized delta function. If 
the approximation holds then the average of the 
conserved current should be closed to its exact 
value. 
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For j3 = 5.1, we get the conserved charge 
to be (Q{t))z c = 3.0004(4) for k = 3 and 
for k = and 6. We will first note that for 
k = 3 we get the conserved charge to be very 
close to the exact value, indicating that there is 
very little mixing with the other allowed sectors 
k = ... — 21, —9, 15, 27.... The other two values are 
due to symmetry. For k = sector, it is easy to 
understand. For the k = 6 sector, the explana- 
tion lies in the fact that N = 12. Thus for every 
allowed sector k + mN there is an m' such that 
k + m'N = — (fc + mN). In that case the sum in 
the numerator vanishes and the conserved charge 
ends up being zero. 

We conclude that there is very little mixing 
with the other allowed sectors in our runs. We 
attribute the smallness of the deviation in the 
case k — 3 to the fact that the free energy for 
our baryons is very large. However, we expect 
that as we lower the quark mass the deviation 
will become more significant. If the deviation is 
too large then we might have to increase the pa- 
rameter N for the discrete Fourier transform. 

8. Conclusions and outlook 

Summing up, we showed that the canonical 
partition function approach can be used to inves- 
tigate the phase structure of QCD. We are able to 
explore temperatures around the critical temper- 
ature and much higher densities than those avail- 
able by other methods. We also find evidence 
that sign fluctuations might limit us in reaching 
much lower temperature and large fermion num- 
ber. We checked that the discrete Fourier ap- 
proximation to the canonical partition function 
introduces only minimal deviations. 

In our runs, we find that even for a small 
lattice, the picture that emerges for the phase 
structure is consistent with expectations. The 
Polyakov loop shows signals of deconfining phase 
transition both in temperature and density direc- 
tions. To complement this picture we would like 
to see a signal in the fermionic sector. For this 
we need to go to smaller masses and, maybe, to 
finer lattices. 

Another goal would be to locate certain points 
on the phase transition line. For this we need to 



get to even lower temperatures and/or densities. 
While, we might be limited in going to much lower 
temperatures, we should be able to reach lower 
densities. 

To reach a lower density with larger volume, 
we need to introduce an esitmator for the deter- 
minant. We should stress that the method used 
to generate the modified ensemble has no bear- 
ing on whether we have a sign problem or not. It 
is an intrinsic property of the modified ensemble. 
Consequently, we do not expect any difference in 
the sign behavior as we introduce the estimator. 

REFERENCES 

1. K.F. Liu, Talk given at 3rd International 
Workshop on Numerical Analysis and Lattice 
QCD, Edinburgh, Scotland, 30 Jun - 4 Jul 
2003, |hep-lat/0312027| 

2. K.F. L iu, Int. Jour. M od. Phys.BW, 2017 
(2002), hep-lat/0202026. 

3. M. Faber, O. Borisenko, S. Mashkevich, and 
G. Zinovjev, Nucl. Phys. B444, 563 (1995). 

4. E. Dagotto, A. Moreo, R. Sugar, and D. Tous- 
saint, Phys. Rev. B41, 811 (1990); A. Hasen- 
fratz and D. Toussaint, Nucl. Phys. B371, 
539 (1992); J. Engels, O. Kaczmarek, F. 
Karsch, and E. Laermann, Nucl. Phys. B558, 
307 (1999), |hep-lat/9 903030; M. Alford, A. 
Kapustin, and F. Wilczek, Phys. Rev. D59, 
054502 (2000). 

5. M. Alford, Nucl. Phys. B (Proc. Suppl.) 73, 
161 (1999). 

6. See for example B. Joo, I. Horvath, K. F. Liu, 
Phys. Rev. D67, 074505 (2003); A. Alexan- 
dra and A. Hasenfratz, Phys. Rev. D66, 
094502 (2002). 

7. M. Faber, O. Borisenko, S. Mashkevich, and 
G. Zinovjev, Nucl. Phys. B(Proc. Suppl.) 42, 
484 (1995). 

8. S. Kratochvila and P. de Forcrand, 
|hep-lat/ 0309146 

9. C. Thron, S.J. Dong, K. F. Liu, and H. 
P. Ying, Phys. Rev. D57, 1642 (1998), 
hep-lat/9707001 



